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Abstract 

The dynamics induced while controUing quantum systems by optimaUy 
shaped laser pulses have often been difficult to understand in detail. A 
method is presented for quantifying the importance of specific sequences 
of quantum transitions involved in the control process. The method is 
based on a "beable" formulation of quantum mechanics due to John Bell 
that rigorously maps the quantum evolution onto an ensemble of stochas- 
tic trajectories over a classical state space. Detailed mechanism identi- 
fication is illustrated with a model 7-level system. A general procedure 
is presented to extract mechanism information directly from closed-loop 
control experiments. Application to simulated experimental data for the 
model system proves robust with up to 25% noise. 

1 Introduction 

Advances in pulse shaping for ultrafast lasers, fast detection techniques, and 
their integration via closed-loop algorithms have made it possible to control 
the dynamics of a variety of quantum systems in the laboratory. Excitation 
may be either in the strong or weak field regime, with the goal of obtaining 
some desired final state. Success in achieving that goal is gauged by a detected 
signal (e.i?., the mass spectrum in the case of selective molecular fragmentation), 
and this information is fed back into a learning algorithm which alters the 
laser pulse shape for the next round of experiments. High duty cycles of ^ 0.1 
seconds or less per control experiment make it possible to iterate this process 
many times and perform efficient experimental searches over a control parameter 
space defining the laser pulse shape. 

As an example of this process, experiments have employed closed-loop meth- 
ods for selective fragmentation and ionization of organic and organometallic 
[|| Q compounds, as well as for enhancing optical response in solid-state and 
other chemical systems Q Yields of targeted species are typically 
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enhanced considerably over those obtained by non-optimized methods. It is 
found that the optimal pulse shapes achieving these enhancements can be quite 
complicated, and understanding their physical significance has proven difficult. 
The same general observations also apply to the many optimal control design 
simulations carried out in recent years Q Q Q Q. 

The present paper will address the identification of control mechanisms in 
theoretical calculations as well as for direct application in the laboratory. In 
we will first describe John Bell's beable model for finite dimensional Hilbert 
spaces, in order to obtain a precise (but non-unique) definition of "mechanism" 
for quantum systems in terms of trajectories over their associated classical state 
spaces. For instance, in molecular systems a trajectory would take the form of 
a sequence of transitions that starts with a given initial molecular configuration 
and switches to another configuration at a distinct time ti , and then to another 
at t2, etc. — to be contrasted with a continuously changing superposition of many 
such configurations. 

The means to numerically implement this mechanism concept is presented 
in An application to the problem of population transfer for a model 7-level 
system is given in which illustrates the usefulness of mechanism information 
in understanding control processes. We then show how the beable approach 
leads to the laboratory working relations (^7|) and ( p^ ) , which make it possible 
to identify some basic aspects of control mechanisms directly from experimental 
data. We illustrate this process in §^ on simulated experimental data for the 
model 7-level system. The overall laboratory algorithm for extracting control 
mechanism information is condensed into a general-purpose procedure in fJ^. 



2 Beables and quantum theory 

Consider a control problem posed in terms of the quantum evolution 

ih:^^m)) = Hit)\m) (1) 

over a finite dimensional Hilbert space with basis \n) where n = 0, 1, 2, . . .. Here 
H{t) = Hq — iJ,E{t) incorporates the effect of the control field E{t), and we can 
explicitly follow the evolution of into a desired final state \4>{t{)). 

This paper is concerned with the question: what is the importance of a 
given sequence rii — > n2 — > • ■ • of actual transitions — or, more specifically, of a 
given trajectory defined as a continuous function n(t) of time — in achieving the 
desired final state |'0(if))? In other words, it is clear that the system is being 
driven into a desired state, but can we find a physical picture of how this is 
being accomplished? 

A conventional answer to the question raised above, essentially that given by 
Bohr on first seeing Feynman's path integral, is to reject the question as ill-posed 
because quantum mechanics is said to forbid consideration of precisely defined 
trajectories over the classical state space {n}. Nevertheless, it is well established 
that there exist dynamical models generating an ensemble of trajectories n(t) 
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whose statistical properties exactly match those associated with \tl}{t)) at each 
t. In the case of a continuous state space, the first such model was that of 
de Broglie, later rediscovered and completed by Bohm [^. They reintroduce 
classical-like particle trajectories into quantum theory by taking the probability 
current J[V'] to describe a statistical ensemble of real particles. So, 



gives the velocity of a particle with mass m and position x at time t, in de 
Broglie-Bohm (dBB) theory. The physical particle is taken to exist indepen- 
dently of, but also to have its motion determined by, the wavefunction ip. The 
time evolution of -0 itself is just given by the Schrodinger equation. 

Bohm developed a full account of how ensembles of such classical-like par- 
ticles could reproduce the predictions of quantum mechanics. A basic issue is 
to compare t) with the statistical distribution P(x, t) describing an ensem- 
ble of particles evolving by (^. One can show that if the initial distribution 
of particles satisfies P(x, 0) = |?/'(x, 0)|^, then P(x, t) — |-0(x, t)|^ will hold 
for all t > 0. That is, if the ensemble is initially in the "quantum equilib- 
rium" distribution given by |'0(x, 0)p, the dynamics — (||) for the particles, and 
the Schrodinger equation for ip — will preserve this equilibrium, consistent with 
the predictions of standard quantum theory. The result is easily generalized 
to arbitrary interacting iV-particle systems by taking x as a point in the SN 
dimensional configuration space. 

In dBB theory the position representation has a special status. While one 
may still regard ■0 as a basis-independent object, the particle dynamics is given 
by (H) specifically in terms of (x|-0) rather than {p\ip) or some other represen- 
tation. But, it is easy to formulate analogs of dBB theory in different bases. 
For instance, one might choose the momentum values p as the beables|^ of the 
theory, and the dBB trajectories x(i) would be replaced by momentum space 
trajectories p{t). 

In the context of a finite dimensional Hilbert space with basis \n), the beables 
can be taken as the sites n of the classical state space {n} analogous to {x} or 
{p}. Some law analogous to (|^) must be given to generate beable trajectories 
n{t) over the state space. Such trajectories would provide a physical picture of 
the quantum transitions induced by a control field E{t). 

In one such theory due to John Bell iQ, trajectories arise from beables 
stochastically jumping between sites connected by non-zero Hamiltonian matrix 
elements. To define this Broglie-Bohm-Bell (BBB) theory, the probability for a 
beable to jump from site m to a distinct site n, sometime in the interval (t, t + e), 
is taken as 



^Bell used the term "beables" rather than the misnomer "hidden variables" to distinguish 
them from conventional observables. 





Tnmit)e = 



2Re{z„„(t)}e if Re{z„™(i)} > 
if Re{z„™(i)} < 



(3) 
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where 

h 1pm[t)* 

and ipn = (^IV'); etc. To ensure normalization, the probabihty for a beable to 
stay at m is thus given by 1 — T„„i(i)e, where the primed sum excludes the 
diagonal term n = m. From (|), we find 

Re{z„„J = -~Re{Zmn} I " 12 (5) 

which implies through that either Tnm{t) — or Trrin (^) — at any partic- 
ular time t. Together with (|l|) this gives 

^IV"™!^ = ^ 2Re{z„„|V'm|^} = y^(J'nm|'0m|^ " T'mn|V'«|^) (6) 
m 7n 

as a type of master equation. Note that the T„m term contributes when 
Rej-^nm} > 0, and the Tmn term contributes when Re{z„m} < 0. 

Now consider the probability distribution Pn{t) of beables in state space 
generated by the jump rule (||). Accounting for the influx and outflux of beables 
at site n, we see Pn{t) satisfies 

m 

which is formally identical to (|6|). Thus, provided P„(0) = |'!/'„(0)p, we are 
guaranteed Pn{t) — |'0n(i)P for all t > 0, which expresses the equivalence of 
BBB theory and ordinary quantum mechanics in terms of statistical predictions. 

The answer to the initial question regarding the importance of a given tra- 
jectory in achieving the desired state |'!/'(if )) is now very simple. The importance 
may be taken as just the probability of realizing that trajectory with the jump 
rule (H) . We can express the final state population in terms of these path prob- 
abilities via the integral (path sum) version of (0): 

PnAU) = 5]P„o(0)Prob(7') (8) 

V 

Prob(p) = n ^^".+1". n (1 - 

where Prob('P) is the probability of realizing the path V = (rio,ni, . . .) under 
(^, and Hp gives the beable configuration at t = pe. The first sum is taken over 
all such paths ending on rif at t = tf, and J[P] = {p \ rtp+i 7^ T^-p} defines the 
jump set. 

The above argument for the equivalence of BBB theory and ordinary quan- 
tum mechanics ensures that the path probabilities Prob(7') are consistent with 
the quantum distribution governing observables. But, it should be noted 
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that BBB theory is not unique in this regard. The rule (^) may be altered in 
non-trivial ways while preserving the master equation (^) The definition 

(||) might even be changed in ways that do not preserve (|), if one is willing to 
relinquish a strict probability interpretation for the trajectories p3| . 

In general, there are many different ways to assign probabilities to trajecto- 
ries that all result in the same time-dependent occupation probabilities P„(t). 
The predictions of quantum mechanics, therefore, cannot select a single assign- 
ment. This non-uniqueness at the root of quantum mechanism identification 
can be dealt with only by reference to the simplicity and explanatory power of 
a given mechanism definition. Below we adopt the definition (M). 



3 Simulating beables in quantum control 

An ultimate goal is to obtain dynamical mechanism information directly from 
experimental data associated with the closed-loop control field optimization, 
without pre-existing knowledge of the system Hamiltonian or wavefunction. 
Methods employing BBB theory for this purpose are presented in §|^, but first 
we shall study control mechanisms for a model system whose Hilbert space and 
quantum evolution are given explicitly in numerical simulations. 

Consider a quantum-optical system with level energies huj„ and dipole mo- 
ments Unm- Applying an external laser field E{t), the Hamiltonian in the inter- 
action picture ||l^ is 

Hi = E{t) J2 Mnme'-^-^ln) (m| (9) 

n m 

where w„„i = "-^n ~ ^m- We will drop the subscript / from now on. E(t) is 
assumed to be given by an independent optimization algorithm designed to, for 
example, maximally transfer population from to \n{). 

A simple second-order Schrodinger propagator was used to solve (|^) in the 
interaction picture, relying on a factorization of the evolution operator as 

N-l , J . 

where tp — pe = pt{ /N and T is the time-ordering symbol. Choosing a time step 
e <C fi./ fJ.E, we can approximate (|l^) by dropping the T operations on the right 
hand side and computing the integrals directly. In doing this an error is accrued 
per time step given by the Baker-Hausdorf identity e"^"*"^ = e'^e^e~2 [^>^1+' ' as 



fj"^ fj^\H{r),H{s)]drds ^ (^)^'^- (H) 



The right hand estimate is obtained by expanding H{r) to first order about 
r = s and noticing that the E'{s) term in H'{s) commutes with H[s). The 
error (pT|) would generally dominate third order terms like (fiEe/h)^ . 
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If the control field is given as E{t) = Re{^j aiEi{t)}, where 



with A{t) and (l){t) possibly adiabatic, we can evaluate / H{s)ds by writing 

(Simply writing J H{s)ds w eH(tp) is not appropriate because we do not want 
to exclude weak field excitation, i.e. fiE <C huj^ so that we ~ 1 may hold.) Thus 
in the adiabatic case |i/'(i)) can be propagated in steps determined by A{t) and 
(f){t) rather than the phase factors e"^*. 

Consider the evolution of beable trajectories according to (||), which appears 
to require a time step small enough that each part of H , including the e'"* terms, 
not vary much over the step. Nevertheless, the total probability of jumping from 
m to n over [tp^tpj^i) is given by the integral J Tnm{s)ds over that range with 
~ {pEe/h)"^ corrections. Thus we can take an effective jump probability for the 
interval {tp,tp^i) as given by (^ with 

Znm{tp) « - ^"^/Z w T- / Hnra(.s)ds (13) 



ipm{tp)* he 

evaluated using (p^). If we <C 1 does not hold, care must be taken to extend the 



integration in (13) only over t e {tp, tp+i) for which Re{z„m(t)} > 0, leading to 
additional boundary terms in the phase difference part of (|l2|). Moving the ^l>* 
ratio outside the integral in ( p^ ) produces an error per time step of order 

e^Hd^J ffiEe 



which is again comparable to (|ll[). Therefore beable trajectories may be prop- 
agated in steps determined by A{t) and (t>(t), i.e. in sink with the Schrodinger 
propagator. 



4 A model 7-level system: mechanism analysis 
of an optimal control design 

The beable trajectory methodology for identification of control mechanisms will 
be illustrated with a 7-level system where u;„ and are given in Fig. |^. The 
(non-adiabatic) control field E{t) shown in Fig. [2| i s obtained from a steepest 
descents algorithm over the space of field histories |14| . It is optimized to transfer 
population from the ground state |0) to the highest excited state |6). By t = 100 
fs, the transfer is found to be completed with approximately 97% efficiency (see 
Fig. I). 



6 



Together with the second-order Schrodinger propagator, using time step e = 
.025 fs, an ensemble of A'traj = 10^ bcable trajectories is evolved, all starting 
in the ground state n = at t ~ 0. At each time step, a given beable at site 
m is randomly made either to jump to a neighboring site n ^ m according to 
the probabilities Tnmi given by (||) with (13), or else stay at m. Four sample 
trajectories are shown in Fig. ^. As a check, one can count the number of 
beables residing on each site n at time t to estimate the occupation probabilities 
Pnit) and verify that they match the quantum prescriptions \ipn{t)\'^- The finite- 
ensemble deviations are observed to be consistent with a (A^traj)'^^^ convergence 
law. 

About 60% of the trajectories generated are found to involve four jumps, 
and of these the trajectories passing through sites n ~ 2,5 arc noticeably more 
probable than those passing through n — 1,4. 6-jump trajectories comprise 
about 30% of the ensemble. And it becomes increasingly less likely to find tra- 
jectories with more and more jumps. The largest number of jumps observed in 
a single trajectory was 14. Three such trajectories occurred out of the ensemble 
total 10^. 

A natural expectation is that the optimal field E{t) would concentrate on 
the higher probability trajectories and not waste much effort on guiding highly 
improbable trajectories, such as the 14-jumpers, to the target state n = 6, as 
the latter have essentially no impact on the control objective (final population 
of the target state). Interestingly, though, the vast majority of even the lowest 
probability trajectories are still guided to n = 6. Apparently, the optimal field 
is able to coordinate its effect on low probability trajectories with that on other 
trajectories at no real detriment to the latter. We shall come back to this point 
later. 

One way to conveniently categorize the large set of trajectories, each express- 
ible as a sequence of time-labeled jumps (ii,ni) (t2,^2) ^ • • is to drop 
the time labels, leaving only the "pathway" iii —> n2 — > • ■ ■. The importance 
of a given pathway is then computed as the frequency of trajectories associated 
with that pathway. Table |l| lists some important and/or interesting pathways 
and their probabilities. 

Fig. H shows some typical trajectories associated with the first and fifth 
pathways listed in Table |^ — involving 4 and 6 jumps respectively. E(t) guides 
the 4-jumpers upward in energy, and they begin to arrive at n = 6 around 
t = 80 fs, early enough that stragglers can catch up but too late for the over- 
achievers of the group to head off elsewhere. This corresponds to the onset of 
heavy growth for \ipQ{t)\'^ around t = 80 fs (see Fig. ||). The 6-jumpers first 
reach n — 6 around t = 50 fs, but almost all fall back to n = 5 by t = 80 fs, 
reuniting with the 4-jumpers just as they begin to jump up to n = 6. These 6- 
jumpers, along with other high-order contributions, thus explain the small surge 
in |'(/;g(t)p between 50 and 80 fs. Another much smaller surge around t = 30 fs 
and one still smaller around t = 20 fs (see inset of Fig. ||) are attributable to 
8-th and higher order trajectories "ringing" back and forth on 5^6. 

For t £ (70 fs, 80 fs), many of the 6-jumpers are at n = 6 and need to 
be de-excitcd on the 6^5 transition before they can jump back up to n = 
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6. Simultaneously, many of the 4-jumpers are at n = 5 and should not be 
prematurely excited on 5 — > 6, lest they not remain at n = 6 through t = 100 
fs. The optimal field thus faces a conundrum: how to stimulate the 2^6 
transition preferentially for the 6-jumpers (in n = 6) over the 4-jumpers (in 
n = 5). The means by which this feat is accomplished may be understood by 
reference to the jump rule (||). E{t) induces jumps through the explicit Hnm{t) 
factor but also through the ip* quotient, which depends on E(t) through ([^). In 
particular, ^ implies that at any one time t jumps on this transition must be 
either all upward or all downward. The active direction is switched back and 
forth according to the sign of Re{z65(t)}. 

Fig. H plots \E{t)\ and Re{z65(i)}, which controls the upward jump rate 
Te5{t). For t £ (70 fs,80 fs) one sees that when \E{t)\ is large, most often 
R-e{z65(i)} dips below zero, disallowing any upward jumps. The correlation 
coefficient between \E{t)\ and Re{ze5{t)} in this range is —0.4955. On the other 
hand, the correlation between \E{t)\ and Re{z5e{t)}, which controls downward 
jumping, is +0.4475 over the same range. 

Looking at the trajectories in more detail, one notices a distinct bunching of 
jumps. Beables tend to jump together in narrow time bands, or else to abstain in 
unison from jumping. This behavior can be gauged by calculating the two-time 
jump-jump correlation function: 

p=0 

where Jo(t) is the number of jumps of type Q occurring in {t, t + e), and is a 
subset of the entire ensemble of trajectories. For instance, the two-time function 
with Q taken as the set of jumps on the 5^6 transition is plotted in Fig. ^ 
The fs time-scale oscillations correspond to the level splittings uJnm and the 
dominant frequency components of E(t). Enhanced correlations around r = 
correspond to the jump bunching noticeable in the trajectories. Two side-bands 
around r = ±40 fs are associated with 6-jump and higher order trajectories 
that go up, down, and up again on 5 <-> 6 over the approximate time window 
(50 fs,90 fs). This conclusion can be verified by computing two-time functions 
with Q specialized to particular pathways. Other much smaller features for 
|t| > 60 fs (see inset of Fig. |^) are attributable to higher order trajectories 
ringing on 5 <-> 6. 

In general, the fs oscillations characteristic of these two-time functions show 
that E(t) works in an essentially discrete way, turning on the flow of beables 
over a given transition and then turning it off with a duty cycle of ~ 2 fs. The 
associated bandwidth of « 0.5 fs~^ is small enough to discriminate between all 
non-degenerate ujnm except between W35(= UJ34) and LiJ5ei= ^45), which differ 
by only 0.12 fs^^. This circumstance leaves effectively three distinguishable 
transitions. With a total time of 100 fs, the control field E{t) can potentially 
enact roughly 150 separate fiow operations. The fact that trajectories with 
pathway probability ^1% are still almost always guided successfully to rt = 6 
suggests that these ~ 150 operations are more than necessary to obtain the 
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97% success rate achieved by the optimal control algorithm in this simulation. 
It appears that the algorithm actively sweeps these aberrant trajectories back 
into the mainstream so as to maximize even their minute contribution to the 
control objective. 

5 Control mechanism identification in the labo- 
ratory 

Using these beable trajectory methods to extract mechanism information di- 
rectly from closed-loop data is complicated by the fact that we cannot assume 
knowledge of a time-dependent wavefunction, Hamiltonian, or possibly even the 
energy level structure of the system. Frequently in the laboratory, the only avail- 
able information consists of final state population measurements and knowledge 
of the control field E{t). 

The following analysis aims to show how a limited statistical characterization 
of beable trajectories may be generated from laboratory data associated with 
a given optimal control field. In particular, we will show how to extract jmin, 
the minimum number of jumps necessary to reach the final state Uf from the 
initial state ni; also {jv), the average number of such jumps over an ensemble 
of beable trajectories; and possibly higher moments {(jv)^) as well. After a 
general formulation of this analysis is presented, it will be applied to simulated 
experimental data in the case of the model 7-level system considered above. 

We propose to obtain mechanism information by examining the effect on 
the final state population IV'nfl^ of variations in the control field away from 
optimality. Consider the simplest such scheme, wherein the amplitude of the 
control field is modulated by a constant M. independent of time: 

E{t) E{t) = ME{t) 

giving rise to a new time-dependent solution \tl){t)) — in particular, a new final 
state population |V-'nf(if)P and new path probabilities Prob(7^). These quanti- 
ties are obtained by taking r„„i — * Tnm in (^, which is to say using E{t) and 
in the jump rule (||). 
To express Prob('P) in terms of Prob(7'), we can write 

n ^"p+i«p = -^^"^ n ^"p+i"p n ^ 

peJ peJ peJ ^ 

T-r \'4'n,+ ^{tp)\ I -pr |^n,+ ifa)| \ 

where j-p is the number of jumps in V and 

7 _ ( .rj V'np+lfa) ^ 

<t)p = arg - . 

V Vnp[tp) J 
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To simplify , note that if j-p were very large, then successive terms in each of 
the last two products would tend to cancel, leaving only endpoint contributions. 
Making the reasonable approximation that they do completely cancel yields 



1^ 

Further, we can make the expansion 



pe.J ^°^^P 

about A4 = 1, where the a^'' depend on the path V but not on A4. And 
similarly: 



Combining these expansions gives a relationship between the path probabilities 
Prob(7') in the modulated case to those, Prob('P), in the unmodulated case, 
which are the ones containing mechanism information regarding the actual op- 
timal control field E{t). We can thus write the final population as 

- 5]|^„„(o)|2p?5b(p) 

V 

where a-p = a^'' + and higher order terms in the expansion have been 
dropped. (This approximation is not as crude as it might seem, since for small 
M away from 1, the behavior of \4'nt{tf)\^ is dominated by the M^'^ factor.) 
Cancelling one power of (tf)|, and recalling that the sum is taken only over 
paths ending on n = nf so that li/^nt (if)P = 'Ylv Prob('P), we have 

\i^nAU)\ « |^n,(if)|(A^^-e-''-(^-i)) (16) 

where (• • •) denotes an average over the trajectory ensemble generated by the 
(unmodulated) optimal field E{t). Beables in this ensemble are taken as initially 
distributed at t = according to |'(/'n(0)|^, and only trajectories that successfully 
reach n = nt at i = tf are counted. 

Note that for ^A close enough to 0, the minimum value jmin taken on by j-p 
will dominate the expectation value in ([l^), and 

log \i^nAil)\ = JminlogA^ + 0(1) (17) 
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gives the dominant behavior independent of a-p. If we suppose that a-p, where 
it is relevant, depends primarily on the endpoints of V, which are fixed, and 
only weakly on the rest of the path, then a-p can be approximated by some 
characteristic value a. Putting Ai^^ = e-'^ log^ in (jj^) and expanding in powers 
of log M now gives 

\Mtf)\ « |^„,(if)|e-''(^-i)f;^^^(logA^)'= (18) 

k=0 

for the final state population under a modulated field, expressed in terms of 
the desired statistical properties of the trajectory ensemble under the optimal 
field itself. Here, a enters as an additional parameter that must be extracted 
from the data. Equations (^7|) and (|l8|) form the working relations to extract 
mechanism information from laboratory data. 



6 Simulated experiments on a 7-level system 



In order to extract quantities like (j-p) using the results (y_7|) and (|8|) data 
must be generated for the final state population IV'nf (if)P at many values of the 
modulation factor M over some range (A^min, Almax) ^ (0,1.5). The desired 
quantities are obtained as parameters in fitting (^7|) and (|8|) to the data as a 
function of Ai. 

One set of simulated data for the above 7-level system is shown in Fig. ^; the 
sampling increment is AM = .01. Noise has been introduced by multiplying the 
exact |'0„j(tf)p values by an independent Gaussian-distributed random number 
for each value of A4, where the distribution is chosen to have mean 1, and 
various standard deviations a have been sampled. 

We can determine jmin from the data using (p7|), which implies 



. .^db^^^ (19) 

For instance. Fig. ^ plots the derivative in (^9|), calculated with finite differences 
from the |'(/'nf(if)P simulated data for a — .1, which correctly gives jmin = 4 
as the limiting value. Determination of jmin proved robust to multiplicative 
Gaussian noise up to the 40% level {a = A). 

The quantity (j-p) is more difficult to extract, because while the sum in 
converges to as — > 0, the terms of the sum individually diverge and must 
cancel in a delicate manner. Therefore truncating the sum to an upper limit 
fcmax becomes a very bad approximation near = 0. This unstable behavior 
can be controlled by carefully setting the range (A^min, A/max) of data to be 
fitted, given a choice of fcmax- 

It is also convenient to constrain the fit by the previous determination of 
jmin — 4. We have done this by noting that if E{t) is truly optimal, then 
|'0rat(^f)l must have a maximum at A4 = 1, which implies that a ~ (jp). This 
can be used as a weaker constraint on the auxiliary parameter a by just requiring 
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a > jmin = 4 in the fit without necessarily supposing that E{t) is exactly 
optimal. We then check that a « (j-p) is satisfied in the fit. Fig. || shows one 
such fit where the fitting range is M £ (.44, .92). One can see that the fit closely 
tracks the data for M in this range but quickly diverges from the data just 
below A4 = .44 (and, less severely, above A4 = .92) due to the sum-truncation 
instability mentioned previously. 

In order to identify appropriate ranges in general, we have searched over all 
combinations such that 

•2 < Mmin < -8 

.7 < A^max < 1.6 (20) 
M max M mill > 10 

Mathematica's implementation of the Levenberg-Marquardt non-linear fitting 
algorithm was used on simulated data for each value of a between and .5 with 
a .01 increment. The best fit at each a was used to determine the value of (j-p) 
most consistent with the simulated data at the given noise level. 

For this analysis fcmax — 4 was chosen somewhat arbitrarily to balance com- 
putational cost and precision. In practice it is likely that the moments {{jv)^) 
for lower k values will be most reliably extracted from the data, especially con- 
sidering the laboratory noise. In the simulations it was found that {j-p) could 
be reliably extracted, but higher moments were unstable and unreliable. For 
example, {{j-pY) was frequently found to lie slightly below the corresponding fit 
values for (j-p)'^, which is inconsistent with the interpretation of these values as 
statistical moments of an underlying random variable jp. Further constraints 
could be introduced to attempt to stabilize the extraction of higher moments, 
but care is needed so as not to overfit the data. 

Fig. |l^ shows the (jp) values obtained by fitting data with a = .1 for each 
choice of (A^min, A^max) and Fig. |l| shows the corresponding quality of each 
fit as measured by its mean squared deviations. In Fig. as well as in the 
corresponding plots for all other values of a studied, two diagonal strips emerge 
running above a set of smaller islands. The surrounding white "sea" comprises 
fits that give (j-p) < 4, which we know to be ruled out by the determination of 

Jmin • 

A virtually identical pattern arises in the fit quality plots. The two strips 
and underlying islands are seen to give much better fits than the white sea. An 
additional connected region of good fits is found to extend across the lower-left 
corner of Fig. nearly all of which are ruled out by jmin — 4. This connected 
region is somewhat pathological because much of it corresponds to fitting ranges 
that fail to capture the important behavior of |-0nt(^f)P near = 1, and 
therefore can be ignored. Then the best fits for all values of cr sampled are 
found to come from the cluster of islands at A4„iax ~ -95. As a is increased 
from to .5, these islands flow from A^min ~ -5 to A^min ~ -3, carrying with 
them the best fit site. Note that the small triangular area in the lower right 
corner, most noticeable in Fig. is a region excluded from consideration by 
the third constraint in (|20|). 
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The best fit values of {j-p) are shown as a function of a in Fig. |T^. These 
values are to be compared with the exact result {j-p) = 4.907 obtained from the 
trajectory ensemble calculations in which require explicit knowledge of the 
level structure and dipole moments /i„m of the system. The ramping behavior 
in Fig. |l2| results from the sampling increment AA^ of the simulated data. 
Transitioning between one ramp and another corresponds to the shifting of the 
best fit location by one or two units of AAl. 

These {j-p) values are in good agreement (3% discrepancy) with the exact 
value for noise at the level of 0-25%. It should be noted that a qualitative change 
occurs in the case of no noise (cr = 0), where the islands all disappear and the 
strips become extended much further on the downward diagonal. Inspecting 
the fits individually indicates that mean squared deviation does not give an 
adequate measure of fit quality in this special case. This anomaly seems due 
to the fact that, in the absence of Gaussian noise from experimental statistics, 
systematic deviations from ( p^ ) associated with the approximation ([l^) become 
important. 



7 Summary procedure for mechanism identifi- 
cation from control experimental data 

In order to extract the basic mechanism information comprising jmin and {jp) 
from quantum control experimental data, the methods of §^ can be distilled into 
the following general procedure: 

1. Perform a closed-loop optimization of population transfer, giving an asso- 
ciated optimal laser pulse shape E{t). 

2. Apply a modulated field E{t) = A4E{t) to the system and measure the 
the resulting final state population \4'ntiti)\'^ for many values of A4 over 
some range ~ (5, 1.5), where (5 is a positive value near determined by 
experimental sensitivity. 

3. Extract jmin from the data by extrapolating the limit (^9|). 

4. Choose a truncation of the sum in (]l8|), e.g. fcmax = 4, and perform a non- 
linear fit to the data for each of a set of fitting ranges (A^min, A^max), £-9- 
the set ( p0| ) . One may choose to constrain the fit by requiring a > j^un in 
dl). 

5. Plot the fit values of {jp), as in Fig. |l^, and the corresponding mean 
squared deviations, as in Fig. ^ over the A^min--A4max plane. Exclude 
regions in which the fit violates the condition {jp) > jmin and pathological 
regions like in the lower-left corner of Fig. |l^. 

6. Find the point (A^min, Almax) at which the mean squared deviation is 
minimized, giving the associated fit value of {jp) as that most consistent 
with the data. 
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8 Remarks and Conclusion 



Since Bell's model can be defined for any choice of basis there is a more 
general question of how mechanism analysis varies with the choice of basis. 
Beyond that, Bell's jump rule (||) itself permits generalization providing 
additional freedom over which trajectory probability assignments may vary. The 
import of this freedom for mechanism identification remains to be determined. 

This paper has shown how Bell's bcablc model of quantum mechanics can 
be used to understand the dynamics of quantum systems driven by complicated 
optimal control fields. Beable trajectories are identified with simple physical 
processes effecting the controlled transfer of population from one state to an- 
other, and aggregations of beable trajectories may be used to compute the im- 
portance of different such processes in the dynamics. In the context of a model 
7-level system, numerical simulations reveal four chief pathways and also a host 
of higher order pathways that are collectively significant on the 40% popula- 
tion level. We have shown how the control field sweeps trajectories into these 
pathways by switching on and off beable flow over specific transitions on a fs 
time-scale. 

Beable trajectory methods were then defined in general for extracting statis- 
tical mechanism information directly from experimental data, without requir- 
ing knowledge of a Hamiltonian or even the level structure of the system under 
study. Application to simulated noisy data for the model system produced the 
correct minimum number of quantum transitions in the control process and the 
average number of such transitions to within 3% at noise levels up to 25%. 
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probability 


pathway 
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13 4 6 


0.018 


2 3 5 6 5 6 


0.005 


2 


0.0007 


023564356 



Table 1: The five most probable pathways, followed by the highest probability path- 
way failing to reach n — 6 ai t = 100 fs, and then the highest probability pathway 
involving a topologically non-trivial cycle in state space. The fractional error in the 
pathway probability P is given roughly by (lO^P)"^''^. 
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Figure 1: The model 7-level system \n) with n = 0, 1, . . . , 6. The transition frequencies 
cOnm in units of fs~^ are shown on the right, and non-zero dipole matrix elements Unm 
m units of 10"^° C ■m are indicated by dotted lines. 
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Figure 2: Electric field E{t) in V/A obtained from an optimization algorithm for 
population transfer from |0) to |6) |l4[. 
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Figure 3: Population |i/'6(t)P as a function of time. Detail for small t is shown in the 
inset. 
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Figure 4: One each of the 4, 6, 8, and 10-jump trajectories generated by the jump 
rule are shown with their sites n plotted against time. For viewing purposes, we 
have displaced them a small amount vertically from each other and tilted the jump 
lines slightly away from the vertical. 
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Figure 5: A sample of 20 trajectories each from the pathways 2 3 5 6 and 2 3 5 6 5 6. 
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Figure 6: The optimal field modulus \E{t)\ (dotted line) and Re{ze5{t)} (full line) in 
fs~^ over the range (70 fs, 80 fs). Their anticorrelation causes beables to be preferen- 
tially selected for the downward transition 6 — > 5 over the upward transition 5^6. 
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Figure 7: Jump correlation function Jq^t) associated with jumps on 5 — » 6, plotted 
against the delay time t for the ensemble of 10* trajectories. Detail for large t is 
shown in the inset. 
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Figure 8: The best fit of (|l^) to the simulated \-ip„f (tf)P data (10% noise) as a function 
of Ai; it occurs over the fitting range A4 G (.44, .92). 
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Figure 9: The derivative is calculated from simulated data with noise level cr = .1; its 
limiting value as logM —^—oo gives jmin- 




Figure 10: Fit values for (j-p) over a set of different fitting ranges (A^min, A^min); 
<7 = .1 here. 
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Figure 11: Fit qualities as measured by the inverse of the mean squared deviations 
between the simulated data and the fit; a = .1 here. The highest fit quality appears 
at (.44, .92). 
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Figure 12: Best fit values for (j-p) as a function of the noise level a, compared to the 
exact value (j-p) = 4.907 (gray line). 
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